Numerical Study of the Cahn-Hilliard Equation in One, Two and 

Three Dimensions 



E. V. L. de Mello 

Instituto de Fisica, 
Universidade Federal Fluminense, 
Niteroi, RJ 24210-340, Brazil 

Otton Teixeira da Silveira Filho 

Instituto de Computagao, 
Universidade Federal Fluminense, 
Niteroi, RJ 24210-340, Brazil 

(Dated: February 2, 2008) 

Abstract 

The Cahn-Hilliard equation is related with a number of interesting physical phenomena like the 
spinodal decomposition, phase separation and phase ordering dynamics. On the other hand this 
equation is very stiff an the difficulty to solve it numerically increases with the dimensionality 
and therefore, there are several published numerical studies in one dimension (ID), dealing with 
different approaches, and much fewer in two dimensions (2D). In three dimensions (3D) there are 
very few publications, usually concentrate in some specific result without the details of the used 
numerical scheme. We present here a stable and fast conservative finite difference scheme to solve 
the Cahn-Hilliard with two improvements: a splitting potential into a implicit and explicit in time 
part and a the use of free boundary conditions. We show that gradient stability is achieved in one, 
two and three dimensions with large time marching steps than normal methods. 

PACS numbers: 64.75. +g, 02.30.Jr, 02.70.Bf, 74.80.-g 
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I. INTRODUCTION 



The theory of the phase-ordering dynamics following a rapid cooling down (or quenching) 
through the critical temperature from a homogeneous or disordered phase into an inhomo- 
geneous or ordered state has been studied for decades jl|, Q, S, Q|. This phenomenon is 
known as spinodal decomposition. Part of the fascination of the field is because the or- 
dering does not occur instantaneously. Instead, a network of domains of the equilibrium 
phases develops, and the typical length scale associated with these domains increases with 
the time. In order words, the length scale of the ordered regions growth with time as the 
different (broken-symmetry) phases compete in order to achieve the equilibrium state. One 
of the leading models devised for the theoretical study of this phenomenon is based on the 
Cahn-Hilliard formulation [sj]. The Cahn-Hilliard (CH) theory was originally proposed to 
model the quenching of binary alloys through the critical temperature |5( but it has subse- 
quently been adopted to model many other physical systems which go through a similar 
phase separation P, 0, 0). 

Recently it has been discovered that high critical temperature superconductors have 
an intrinsic inhomogeneous phase which exhibit patternings which involve nanoscale re- 
gions of phase separations, often referred as stripes as revealed by neutron diffraction 
studies j(| and EXAFsQ. Similar findings were measured by very fine scanning tunnel- 
ing microscopy/spectroscopy (STM/S) dataj^, Q| which have shown the local charge and the 
superconducting gap spatial variation through the differential conductance. The ubiquity 
of such intrinsic inhomogeneities, as well as the importance of the materials that exhibit 
them has motivated intense experimental and theoretical research into the details of the 
p_a0flyHQ. Rocont.y a theory of the eritieal Md H a based on a distri- 



bution of superconducting regions with different critical temperatures due to these intrin- 
sic inho mog eneities has explained some nonconventional features of current data on some 
cuprates[15]. The manganites which exhibit a colossal magnetoresistance have also some 
properties which may be linked with clusters formation upon cooling or an intrinsic phase 
separation jl^l. Il7 1 . Researches on these materials are currently been investigated by a sizable 
fraction the condensed matter community. Therefore, studies based on the CH equation may 
be useful to understand the puzzle of the origins of such intrinsic phase separation in these 
materials and how it affects their physical properties. Here we want to study specifically 
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the Cahn-Hilliard equation and we will deal with the problem of phase separation in high-T c 
superconductors in another publication. 

In an alloy system composed by a binary mixture we can define the local phase variable 
u(x) as the difference in concentration between the two (incompressible) components or 
simply, the concentration of one of the components over the domain Q (x G Q) which 
represents the system. It is clear that this type of phase variable or order parameter is always 
conserved for an isolated system and we will use below this property as the main guide for 
our numerical method. In the CH theory, the time variation of the order parameter u(x, t) is 
given in terms of the functional derivative of a time- dependent free-energy functional F(u) 
leading to an equation of motion to diffusive transport of the order parameter, namely 

= MV 2 / 1 

where / is the free-energy density and the (constant) mobility M will hereafter be absorbed 
into the time scale although there are cases in which the mobility can be a function of the 
position jisl]. The free-energy functional is assumed, to most of the physical applications, to 
follow the Ginzburg-Landau (GL) form: 

/ = \e 2 \Vu\ 2 + V(u) (2) 

where the potential V(u) may have a double- well structure like, for instance, V(u) = (u 2 — 

nnnn 

l) 2 /4 as many authors have usedjjlQl, |20, |2l|, |22| and, for the sake of comparison with 
these previous works, we will adopt it below. There are other possibilities like V(u) = 
au 2 /2 + bu A /4 + ... which is more convenient for physical applications since, usually the 
GL free energy is a power expansion in u with coefficients (a, b, ...) which depend on 
the temperature, the applied field, and on other physical properties. It is easy to see that 
the above double well potential will favor two phases with densities u — ±1. If one uses a 
potential with three minima, it will appear three major phases and so on. Bray 4] pointed out 
that one can explore the fact that the order parameter is conserved and the CH equation can 
be written in the form of a continuity equation, d t u = —V.J, with the current J = V(Sf /Su). 
Therefore we may write the CH equation as following, 

^ = _VVV 2 M + M - M 3 ) (3) 
ot 

Normally, < e <§; 1 because the fourth derivative term requires a large stencil and it is 
related with the size of the interface between regions of two different phase. To accurately 



resolve these interfaces a fine space discretization Ax is necessary. The linear term is re- 
sponsible for the interesting dynamics including the instability of constant solutions near 
u = and the nonlinear term is the one which mainly stabilizes the flow[19]. As it has 
already been pointed out |igL I20 , 21, 2^|, both the V 4 and the nonlinear term make the CH 
equation very stiff and it is difficult to solve it numerically. The nonlinear term in principle, 
forbids the use of common Fast Fourier Transform (FFT) methods and brings the additional 
problem that the usual stability analysis like von Newmann criteria cannot be used. These 
difficulties make most of the finite difference schemes to use time steps of many order of 
magnitude smaller than Ax and consequently, it is numerical expensive to reach the time 
scales where the interesting dynamics occur. This is the reason why numerical simulations 
based on Runge-Kutta schemes had to be performed in large supercomputers [2c 

In order to deal with such restrictions, Eyre developed a semi-implicit method [19j] which 
resolves the problems associated with both the stiffness and solvability. Furthermore Eyre 
proved that his algorithms are unconditionally gradient stable for the CH and also for the 
Allen-Cahn equation, which means that the free energy does not increase with time. Both 
gradient stability and the conservation of mass provide us with a simple and rational form 
to establish the stability criterion for the CH equation that replaced the von Neumann 



stability criteria. 
Carlo simulations 



urthermore, mostly either finite difference calculations 
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2J or Monte 



25j | uses periodic boundary conditions as a tentative to mimic a large 



system and we show here that, as concerns phase separation, that free boundary conditions 
are less stiff and more faster achieved. 

The goal of this paper is to make a combinations of a systematic study of the CH equation 
in ID, 2D and 3D using a simplification of the Eyre's method and free boundary conditions. 
To give a better perspective to the approach, we make also calculations with a Crank- 
Nicholson like (CN) implicitly scheme 26], which is unconditionally convergent in ID and, 



using the concept of alternating direction interaction (ADI) method [26j, |2j|, in 2D. These 
calculations demonstrate the advantage of the Eyre's method over the CN scheme. Therefore 
we apply Eyre's method in 3D and study the phenomenon of spinodal decomposition in three 
dimensions. The application to the high T c superconductors and manganites with the study 
of the relevant parameters to their phase diagrams is under current investigation and will 
be discussed in a future work. 
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II. THE PROPERTIES OF THE CAHN-HILLIARD EQUATION 



In order to solve numerically Eq.© by means of a finite difference scheme, we need an 
initial condition over the entire domain, u(x, 0), which usually is a random function or some 
small fluctuations over a specific average and also some type of boundary conditions (BCs). 
During our simulations we have found that these initial conditions, the BCs and the size of 
the system greatly influence the solution u(x, t). The general convenient and flux-conserving 
boundary conditions are^ 



Vu.n\ s& dvt = (V 3 ii).n|xedn = 



(4) 



where n is the outward normal vector on the boundary of the domain Q which we represent 
by dfl. These two equations together are equivalent to 



V(Sf/Su) Sedn = 0. 



(5) 



These BCs lead to the two very important properties which will be our guide to know 
whether our CH numerical solutions are convergent, namely, the conservation of the total 
mass of the system M t ; 

d f f du(x,t) 
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This last equation shows that appropriate solutions of the CH equation must dissipate 
energy and this is called gradient flow. Therefore a time stepping finite difference scheme is 
defined to be gradient stable only if the free energy F(u) does not increase with the time, 
i.e., obeys Eq. ([7j) . Since it is not convenient to use the von Neumann stability analysis, 
gradient stability is regarded as the best stability criterion for finite difference numerical 



solutions of the CH equation 



22j . Furthermore, unconditional gradient stability means that 



the conditions for gradient stability is satisfied for any size of time step and this will be 
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our guide through the simulations. This way to examine the stability of a central difference 
scheme by following the energy was already proposed to nonlinear problems a long time ago 
by Park[28]. We should mention that a similar analysis, based also on the Eyre's approach, 
performing numerical tests of stability and with a very complete classification scheme for the 
stable values of At for the CH and Allen-Cahn equation in 2D was recently developed j^]. 



III. THE DISCRETIZATION METHOD 

Eyre has proposed a semi-implicit method that is unconditional gradient stable if the V(u) 
is the usual two minima potential used in a typical GL free energy and can be divided in 
two parts: V(u) = V c (u) + V e (u) where V c is called contractive and V e is called expansive ^19]. 
He showed that it is possible to achieve unconditional gradient stability if one treats the 
contractive part implicitly and the expansive part explicitly. In our case, since we are using 
V(u) = (u 2 - l) 2 /4, we have: 

u 4 + 1 u 2 
V c {u) = — - — and V e (u) = — — . 

To implement Eyre's scheme we define U^- k {i^j^k = 1,2, ...,N; n = 0,1,2,...) to be the 
approximation to u(x,t) at location x = ih, y = jh, z = kh and t = nK, where h = Ax = 
Ay = Az, K = At and N = Lj Ax. L is the linear size of the system, assuming to be cubic, 
for simplicity. With these definitions, we can write the method for the CH equation as 

Tjn+l _ jjn 

l]k K l3k = -e 2 ^ 1 + V 2 ({U^f - V» h ) . (8) 

In this equation the standard centered difference approximation of the 3D Laplacian operator 
V 2 is 

^ 2 Uijk — iPi+ijk + UJj+lk + Uijk+\ ~ ^ijk + ^i-ljk + Uij-Xk + Uijk-l) I ^ (^) 

which is second order in the spatial step h. The Eq.(JSJ) represents a large coupled set 
of nonlinear equations due to the cubic term. The way to go around this problem is to 
splitting or to linearize it at every time step. Consequently the term (U^ 1 ) 3 is transformed 
into {UJjk) 2 ^t lading the Eq.(jHJ) into a set of linear equations in the step of time n+ 1. It 
has been argued that this nonlinear splitting has the smallest local truncation error [lfll I22I] 
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and therefore it will be the scheme adopted in this work. Then we finally obtain the proposed 
finite difference scheme for the CH equation which is linear (in the above sense), namely, 

Tjn+l _ Tjn 

l]k K l ° k = -e 2 ^ 1 + V 2 ({U^fU^ - U? jk ) . (10) 
or, separating in different times to see the semi-implicit character of the approach, 

+ Kie^U^ 1 + V 2 {U- k fU^) = U? jk - KV 2 U^ k . (11) 



As noted by Furihata et al|20j], the discrete associated boundary conditions equivalent to 
Eq.(j3J), second order in the spatial variable, become 

Tjn —TJ n TT n —TT n 

U ljk — u 5jk J u 2jk — u 4jk \ LZ ) 
^Njk = ^N-Ajk i ^N-ljk = ^N-3jk (13) 

and similar equations for the boundary conditions over the second and third indices repre- 
senting the y and z-directions, respectively. Notice that, differently than most numerical 
simulations applied to physical systems |29l|. these boundary conditions are not periodical. 
Imposing periodical BCs which is common in physical applications and which is largely used 
to artificially simulate a larger system, would bring additional constraint to the solutions 
and, as we show below, the solutions will be more stiff. Therefore the above boundary 
conditions minimize the finite size effects and should be preferably used as shown below. 

The Eq. (jll)) with the above boundary conditions define the finite difference scheme which 
we use in this present work. In the following section, we analyze the numerical results and 
compare with different CN semi-implicit approaches for several dimensions with the same 
double-well potential. 

IV. THE RESULTS OF THE SIMULATIONS 

We start with the study of the CH equation in ID because it is faster then 2D and 3D 
and there are many results which can be used to compare with our simulations. The usual 
semi-implicit and explicit Euler's schemes for the ID CH equation are not gradient stable 
and require very short time intervals At = K, while the implicit Euler's scheme is gradient 



Q0,Q 



stable with K < h 2 /A = 6.25 x 10" 6 [1J, |2l|, |2fl|. The CN-like scheme that we will use 
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for comparison also suffers from this solvability restriction, and requires a minimum time 
interval K < h 2 /9. 

We performed calculations with linear chains of N = 54, 104 and 504 sites and with 
Ax = h = 1/(N — 4), and < x < 1. For all these cases we used K — h/2 which is clearly 
several order of magnitude larger than typical Euler's schemes and, despite this large time 
step, the gradient stability is observed at all times. In Fig.l we show the results for the 
total mass and the free energy F in arbitrary units as function of the running time. In fact, 
using shorter times as shown in Fig.l, we observe that there is an initial transient period in 
the time evolution before gradient stability is achieved. This can be easily seen at at most 
small values of At in Fig.l for either Eyre's or CN-like methods. This transient before the 
stability takes place is connected to the finite size of the system, the initial conditions and 
the BCs. 

As already mentioned, in order to compare our calculations based on Eyre's method with 
other schemes, we have also performed similar simulations with a widely used methodil9u2ll. 

nn 

a semi implicit CN-like scheme 26, 30]. Briefly, the CN method consists in an average of the 
right hand side of Eq.fjHl) at different times: 1/2 at n and 1/2 at n+ 1 what clearly results in 
a semi-implicit scheme in time. According to the above formula, for h = 1/50, this scheme 
will converge for K < 1/(2500 x 9) « 1/20000). Indeed we can see in Fig.l that the CN 
simulations agree very well with the the Eyre's results for the case with K = 1/50000 and a 
similar transient period before the stability is achieved is observed. To compare with other 
studies which use the same parameters for the CH equation, we used an initial condition for 
the linear chain similar to one used by Furihata et aljsjj], namely, 

C/P = 0.1sm(2.7r(i - 3)h) + 0.01cos(47r(i - 3)h) 

+0.06sm(47r(i - 3)h) + 0.02cos(107r(i - 3)h) (14) 

for 3 < i < N — 2 and Uf = 0.03 otherwise (at the borders). This function is shown in 
Fig.©. 

Now that the parameters which gradient stability is established in our system, we study 
the process of spinodal decomposition in ID starting with the initial condition given by 
Eq. (ll4j) . Using the long time step At = 1/100, we see for the different [/" time profile in 
Fig. 2, that at very few time evolution (n = 5) the system shows a tendency to decompose 
in one high and other low density phases. At n — 10 the difference in densities increases 
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1 D Mass-Energy test (a=b=1 , t =0.001, Ax=1/50) 
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FIG. 1: Time evolution of the total mass and the total energy in arbitrary units up to 10 5 time 
steps (abscissa) by Eyre's and Crank-Nicholson's scheme at different time steps for comparison. 

and the system separates in regions which almost reach the limit ±1 values. At n = 50 and 
n = 100 the system is almost entirely separated in the two ±1 phases. At n — 500 the limit 
configuration is reached and the spinodal decomposition is total with the low density phase 
at left and the high phase at right as seen in Fig.2. Furihata et alQ have found a similar 
result starting from the same initial condition using a different method. Larger systems 
behave in a similar fashion with the same sort of decomposition seen in Fig.(J2J). Notice that 
the system separates in two regions with different density (-1 at left near x = and +1 at the 
other and at x = 1). If we had imposed periodic boundary conditions, it would impose an 
additional constraint and the solutions would be different than that two regions of Fig. @ , 
namely, three phase regions with one at the center and two of the same kind at the borders. 
This more complex final configuration takes more time steps to be realized as we checked, 
performing also calculations with periodic boundary conditions. Thus, we demonstrate that 
the used "free" boundary conditions are more "natural" and faster than the periodic ones 
and this is an important finding which will be used in the 2D and 3D studies. 

Now, let's use what we learned above and turn our attention to the 2D problem. We 
worked with systems of sizes of 104 x 104 and 504 x 504 and again Ax = h = 1/(A — 4) 
with N = 104 or 504. We adop t the Eyre's scheme in the form of an alternating dimension 
implicit (ADI) problem 0, 13; that i s > we use d a half marching time step in one direction, 
say, in the x-direction and another half time step in the ^/-direction. Since most finite 
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U(t) (a=1, b=1, £ =0.001, At=1/100Ax=1/50) 




FIG. 2: The time evolution of the profile of U™ for different times starting at the initial condition 
(t = 0) till the total spinodal decomposition is attained at n = 500. The size of the system is 
L = 1. 

difference methods for partial differential equations like the diffusion or heat equation in 2D 



uses the Crank-Nicholson scheme! 26 



27| | , we again made simu 



26 



27 



ations with this scheme which 



is appropriate to be used in connection with the ADI 

In Fig. 3 we plot the results for the 2D mass and energy in arbitrary units. We used 
£ = 0.01 which is slightly bigger than the ID value in order to have larger phase domains 
and we kept the other parameters equal. The initial conditions were chosen to be small 
variation around the average value = 0.5 although variation around zero give the same 
type of result. Thus the initial condition is 



C/P. = 0.5 + Esin(ir(i - 3)20h)sin(n(j - 3)20/i) 



(15) 



for i, j toward the middle and with = 0.5 toward the boundaries of our square systems. 



The choice of a different initial condition than Uf- = is to brake the symmetry and, in 



this way, better identify the formation of the two ±1 density phases (see Fig.(j3J) below). 
Studying the stability conditions, we see that the 2D system requires smaller values of At 
than the ID system. This is because the boundaries are much larger in 2D and the possibility 
to mass and energy flow is greatly enhanced and the instability during the initial transient 
period is larger. For instance, the initial instabilities in the 504 x 504 is even bigger than in 
the 104 x 104 system. 
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The CN scheme requires at least At = Ax/1000 which is around four times the minimum 
ID value and the Eyre's requires at least At = Ax/50. In the simulations with the CN 
scheme, the mass oscillates wildly up to 10 5 steps and the energy has also a small increase 
near this number of marching time steps. The Eyre's results for At = Ax/ 1000 oscillate 
in the beginning transient of the simulations but become gradient stable after 50000 steps 
when the mass and the energy stabilize. On the other hand the simulations are gradient 
stable and does not display any transient oscillation since the beginning for At = Ax/100 
but there is a large loss of average mass, as shown in Fig.®. 



2D Mass-Energy test (a=b=1, e=0.01, Ax=1/100) 
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FIG. 3: The stability of the mass and the energy in 2D for the 104 x 104 system by the Crank- 
Nicholson and Eyre's scheme for Ax = 1/100 and different time marching steps At as specified in 
the legends. 

As mentioned above the value of e must be higher than that used in ID in order to 
observe the phase separation phenomenon, otherwise the phase domains become very small. 
In order to study the phase separation through the 2D Eyre's method (with At = Ax/1000), 
we started with the initial condition of Eq. (jl5|) . As expected from the above above analysis, 
around the 1000 th time step the system start to separate into the ±1 density phases. At the 
1000 t?i time step the be ginning of the spinodal separation is clear as shown in Fig.(j3J) below. 

Around the 6000 t/l time step the the size of the low density islands reach an equilibrium 
size and the spinodal separation is complete but, the boundaries effects are large. We see 
a concentration of the low density phase at the borders and this effect prevents the system 
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FIG. 4: At the top panel, the beginning (n = 2000) of the spinodal decomposition with the 
formation of islands with the values of -1 density in white in a background which tends to +1 
density (in black) because the initial conditions with broken symmetry. At the panel below, we see 
the system at a later time (n = 50000) when it is already in equilibrium. The initial configuration 
is that given by Ea. (fT5J) . 
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FIG. 5: The same type of figure of Fig.Q but for V(u) coefficients a = 0.6 2 and b = 4.0 2 
(instead of a = 1 and 6 = 1, which produces two easily separated phases with u = ±0.15) after the 
equilibrium is achieved (n = 40000). 

to reach a total phase separation as we have seen in ID (see Fig. (2)). This shows how the 
boundary effects are very strong and more significant in 2D. After this time step the system 
achieves a stable configuration without appreciable changes up to 10 6 time steps. Below 
we show a snap shot of this equilibrium configuration at n — 50000. Similar results were 
obtained with the 504 x 504 system but in this case there is a proportional increase of the 
islands size and the equilibrium situation is similar to five times blow up of Fig.(J3J). Simula- 
tions with a larger value of the non-linear coefficient b (see discussion after Eq.(J2}) enhances 
the mass flow inside the system and it is possible to achieve complete phase separation at 
an earlier time, exactly as seen in ID, but with two phases with smaller values than ±1. In 
Fig.© we plot a situation with a = 0.6 and b = 4.0. In this case there is an easy mass flow 
through the system and an state of complete phase separation is reached around n = 40000. 
Notice that the system reaches the equilibrium with nonperiodical BCs. 

We turn now to analyze the 3D system. As we have already found in the 2D case, the 
problem of mass/energy flow is largely enhanced through the boundaries. To deal with such 
effects, we have to use a small (compared with those used in the 2D case) marching time 
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step of At = Ax/10000 with Ax = 1/100 for the (104) 3 system to minimize errors in the 
derivatives near the boundaries. Notice that this time step is much smaller than the one in 
ID and one-two orders of magnitude smaller than in 2D. This is small but still feasible to 
reach the interesting dynamics even in a typical PC of 1GHz. Since Eyre's method is very 
efficient and faster than the majority of other methods used in 3D, this is the only approach 
that we used in 3D. The generalization of basic Euler or Runge-Kutta methods to the CH 
equation in 3D are much more slower, and we can see why there are very few studies of the 
CH equation in 3D. Indeed the values of At ~ 10~ 6 used in our 3D simulations are of the 
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used in ID simulation as we 



same order of magnitude of typical Euler's method 
discussed in the beginning of this section. The adopted initial condition for the 3D system 
is: 



U° jk = 0.01 + esin(ir(i - 3)20h)sin(n(j - 3)20h)sin(Tt(k - 3)20h). 



(16) 



We used an average small initial mass of 0.01 just to minimize the mass flow through the 
boundaries and indeed the mass remains very stable for up to 4 x 10 4 time steps, when it 
starts a slightly increase. 

3D Mass-Energy test (a=b=1 , e=0.03, Ax=1/100) 




20000 



40000 



60000 



80000 



FIG. 6: The stability of the mass and the energy in 3D for the 104 3 system by the Eyre's scheme 
for Ax,y,z = 1/100. 



As concerns the energy, it remains stable for the first step of the simulations but the 
boundaries effects are manifested near the time step 4000 and increases up to step 13000. 
During this transient period of time, although the mass is stable, the simulations are not, 
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as one can conclude by the energy increase in this interval as shown in Fig.©. After this 
time interval (n ~ 13000) the system stabilizes itself and remains gradient stable up to the 
rest of the calculations (up to n — 80000) as can be seen by the decrease of the energy in 
Fig.©. Comparison with the ID and 2D systems reveals clearly that the stability is more 
difficult to be established as the dimensionality increases, as it is natural when one uses any 
finite difference scheme. As Fig.® shows, the system relax from the initial conditions and 
becomes uniform around n = 3000. Around n = 6000 a spinodal decomposition is happening 
and the oscillation in the densities forms an almost uniform pattern as shown in the top 
panel of Fig.©. This figure shows the configuration of the middle plane in the z-direction 
(k = 52) of the 104 x 104 x 104. At the down panel of Fig.© we show a snap-shot of 
n = 8000 which shows the beginning of of the phase separation process as the time flows. 
Above n = 20000 the two phases start to segregate and this segregation process is smooth 
and unconditionally gradient stable as seen in Fig.©. 

At much later time like n = 40000 the phase separation phenomenon is on the way. At 
n = 80000 the phase separation is almost total in the plane passing by the middle of the 
system, as it is shown in the down panel of Fig.©. 
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FIG. 7: After relaxing from the initial condition (Eq. ljlfij) ). we show at the top panel the beginning 
(n = 6000) of the spinodal decomposition with the formation of small islands with the ±1 values 
because the initial condition with broken symmetry. At the down panel, the system at a later time 
(n = 8000) during the process of phase separation. 
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FIG. 8: The continuation of Fig.Q for many time steps later. On the top panel we plot the case 
for n = 40000 and on the down panel the n = 80000. 

V. CONCLUSION 

We have shown in this work that the semi-implicit method due to Eyre, which divides the 
potential V(u) in two parts and treats the contractive part implicitly and the expansive ex- 
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plicitly combined with free boundary conditions achieves unconditional gradient stability in 
ID, 2D and 3D. The method also allows the use of very long marching time steps, compared 
with usual explicit or implicit Euler's schemes, what is convenient because it captures very 
easily the short and long dynamics, characteristic of the Cahn-Hilliard equation. Our sim- 
ulations have demonstrated that gradient stability and spinodal decomposition is achieved 
faster than the normal Euler and Crank-Nicholson methods what is highly deseirable spe- 
cially in 3D where normal methods fail to capture the long dynamics due to the required very 
short time steps. The use of "nonflow" free BCs are more appropriated and also converges 
faster to a final equilibrium configuration than the widely used periodic BCs. 

We believe that the present systematic study may be pertinent to the several branchs of 
physics. The fast scheme developed here are suitable to the study of the spinodal dynamics 
and also to high correlated electron systems with phase separation. We expect therefore, to 
perform these works in the near future. 
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